%end_day = ceil(time_interval/time_step)
%price_path is the matrix where
    %each row is a price path
    %column(i) is the (i*time_step) day
    %max comlumn = end_day+1 because we simulate price_path from t = 0
    %max row = num_of_simulate*2 because we use antithetic variables method
%delta_t = time_step/365 years
%The formula we use to simulate here
%S(t+delta_t) = S(t) + delta_t * f(t,t+delta_t) *S(t) + sigma(t,S_t) * S_t * sqrt{delta_t} * Z
%%
function [discount_factor,price_path] = simulatePricePath(num_of_simulate,time_step,time_interval,initial_price,input_file) 
    %Load data from file
    [maturity_array, zero_coupon_rate, time, fractionS0, volatility] = getData(input_file);
    
    end_day = ceil(time_interval/time_step);
    %Simulate price path
    [discount_factor estimate_risk] = estimateRiskFree(maturity_array,zero_coupon_rate,time_step,end_day);
    
    price_path = zeros(num_of_simulate*2,end_day+1);
    price_path(:,1) = initial_price;

    delta = time_step/365;
    %Z is a matrix of iid ~ N(0,1)
    Z = randn(num_of_simulate*2,end_day);
    %Use antithetic method, we divide Z from 2 part
	Z(num_of_simulate+1:end,:) = Z(num_of_simulate+1:end,:).*(-1);
    
    %Meshgrid to interpolate sigma
    St = fractionS0 *initial_price;
    [St,time] = meshgrid(St,time);
    
    for j=2:end_day+1
        sigma = (estimateSigma(St,time,volatility,price_path(:,j-1),(j-1)*time_step))';
        %Simulate price paths
        %price_path(:,j) = price_path(:,j-1).*(1 + estimate_risk(j-1) + sqrt(delta)*sigma.*Z(:,j-1)')';
        price_path(:,j)=price_path(:,j-1).*exp(estimate_risk(j-1) - 0.5.*sigma.*sigma.*delta + sigma.*sqrt(delta).*Z(:,j-1));
    end
   